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In this work we consider the phase transition from ordered to disordered states that occur in the Vicsek 
model of self-propelled particles. This model was proposed to describe the emergence of collective order in 
0^ \ swarming systems. When noise is added to the motion of the particles, the onset of collective order occurs 

through a dynamical phase transition. Based on their numerical results, Vicsek and his colleagues originally 
concluded that this phase transition was of second order (continuous). However, recent numerical evidence 
seems to indicate that the phase transition might be of first order (discontinuous), thus challenging Vicsek's 
original results. In this work we review the evidence supporting both aspects of this debate. We also show new 
numerical results indicating that the apparent discontinuity of the phase transition may in fact be a numerical 
artifact produced by the artificial periodicity of the boundary conditions. 
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I. INTRODUCTION 

The collective motion of systems such as schools of fish, swarms of insects, or flocks of birds, in which hundreds or thousands 
of organisms move together in the same direction without the apparent guidance of a leader, is one of the most spectacular 
examples of large-scale organization observed in nature. From the physical point of view, there has been a drive to try to 
determine and understand the general principles that govern the emergence of collective order in these systems, in which the 
interactions between individuals are presumably of short-range. During the last 15 years, several models have been proposed to 
account for the large-scale properties of flocks and swarms. However, in spite of many efforts, the understanding of these general 
principles has remained elusive, partly due to the fact that the systems under study are not in equilibrium. Hence, the standard 
theorems and techniques of statistical mechanics that work well to explain the emergence of long range order in equilibrium 
systems, cannot be applied to understand the occurrence of apparently similar phenomena in large groups of moving organisms, 
or systems of self-propelled particles (SPP),as they have come to be known. IHOIl 

One of the simplest models conceived to describe swarming systems was proposed by Vicsek and his collaborators in 1995, 
[Q]]. This model is "simple" in the sense that the interaction rules between the particles are quite easy to formulate and implement 
in computer simulations. However, it has been extremely difficult to characterize its dynamical properties either analytically or 
numerically. The central idea in Vicsek's model is that at every time, each particle moves in the average direction of motion 
of its neighbors, plus some noise. In other words, each particle follows its neighbors. But the particles make "mistakes" when 
■ evaluating the direction of motion of their neighbors, and those mistakes are the noise. There might be several sources of noise. 
For instance, the environment can be blurred, and, consequently, the particles do not see very well the direction in which their 
neighbors are moving. We will call this type of noise extrinsic noise because it has to do with uncertainties in the particle-particle 
"communication" mechanism. On the other hand, it is possible to imagine the particles as having some kind of "free will", so 
that even when they know exactly the direction of their neighbors, they may "decide" to move in a different direction. We 
will call this last type of noise the intrinsic noise because it can be thought of as being related to uncertainties in the internal 
. ,_| | decision-making mechanism of the particles. 

In their original work, Vicsek et al. used intrinsic noise and showed that, when the noise amplitude is increased, the system 
undergoes a phase transition from an ordered state where all the particles move in the same direction, to a disordered state where 
the particles move in random uncorrelated directions (see Fig|TJ. Based on numerical results, this phase transition was thought 
to be of second order because it seemed to exhibit several of the characteristic properties of critical phase transitions, such as a 
non-trivial scaling behavior with the system size, divergence of the spatial correlation length, etc. However, in 2004 Gregoire 
and Chate challenged these results, [2]. By numerically analyzing large systems with many more particles than the ones studied 
by Vicsek and his colleagues, Gregoire and Chate concluded that the phase transition is in fact of first order for both types of 
noise, intrinsic and extrinsic. This elicited a debate about the nature of the phase transition in the Vicsek model, namely, whether 
this phase transition is first order or second order. The existence of this debate, into which some other research groups have got 
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FIG. 1: Snapshots of the Vicsek model, (a) For very small noise (tj = 0.1 ) all the particles move in the same direction, (b) For very large 
noise (77 = 0.8) the particles move in random uncorrelated directions. The shaded circle at the center of each frame indicates the size of the 
interaction vicinity. In both cases intrinsic noise was used. The other parameters are: Number of particles, N = 2000; size of the box, L = 32; 
particle speed, vo = 0.05; interaction radius, ro = 1. 



involved, is indicative of the difficulties posed by the system, for not even such a fundamental property as the nature of the phase 
transition occurring in one of the simplest models of flocking, has been conclusively resolved. 

The question as to whether the phase transition is first order or second order is not just of academic interest. Consider for 
instance a school of fish. For low values of the noise, the nature of the phase transition is irrelevant because the fish peacefully 
move in the same direction in a dynamical state far away from the phase transition. However, when the system is perturbed, for 
instance by a predator, the entire school of fish collectively responds to the attack: The school splits apart into smaller groups, 
which move in different directions and merge again later, confusing the predator. A possible hypothesis for this collective 
behavior to be possible is that, when attacked, the fish adjust their internal parameters to put the whole group close to the phase 
transition. If this phase transition is of second order, then the spatial correlation length diverges making it possible for the entire 
school to respond collectively. However, this collective response would be very difficult to mount if the phase transition was first 
order, since in such a case the spatial correlation length typically remains finite. 

In this work we review the evidence supporting the two sides to this debate, i.e. the evidence supporting the continuous 
character of the phase transition and the evidence supporting its discontinuous character. It is important to mention that this is 
not a review on the general topic of animal flocking. For such a review, we refer the reader to Refs. 

Here we focus only on the debate about the nature of the phase transition in the Vicsek model with intrinsic noise. [41] In the 
next section we introduce the mathematical formulation of this model. In Sec. [Ill] we present the work supporting the conclusion 
that the phase transition is continuous, and in Sec.|IV]the evidence supporting the opposite conclusion. In Sec. [V] we present new 
evidence indicating that the discontinuous phase transition that has been reported may be an artifact of the boundary conditions, 
which in the numerical simulations performed so far have always been fully periodic. Finally, in Sec. [VI] we summarize the 
different aspects of this debate. 



II. THE VICSEK MODEL 



The model consists of N particles placed on a two dimensional square box of sides L with periodic boundary conditions. 
Each particle is characterized by its position x„(t), and its velocity v„(t) = voe ,e "('\ which we represent here as a complex 
number to emphasize the velocity angle 0„(t). All the particles move with the same speed v'o. What changes from one particle's 
velocity to another is the direction of motion given by the angle 9„(t). In order to state mathematically the interaction rule 
between the particles, let us define U„(ro;t) as the circular neighborhood or radius ro centered at x n {t), and as k n (t) the number 
of particles withing this neighborhood (including the n-th particle itself). Let V„(t) be the average velocity of the particles within 
the neighborhood U„(ro;t): 
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FIG. 2: Time series of the order parameter \jf^ (t). Starting from a random initial condition for which yrj (0) ~ 0, after a transient time (about 
2000 time steps in this case), the system reaches a steady state reflected by stationarity of the order parameter. The simulation was carried out 
using intrinsic noise for a system with N = 2 x 10 4 , tj =0.1, voAt = 0.05 and ro = 0.4. 



With these definitions, the interaction between the particles with intrinsic noise is given by the simultaneous updating of all 
the velocity angles and all the particle positions according to the rule 



9„(t+At) = Angle |v n (f)J +T}§ B (r), 

v„(f + Ar) = v oe '- e " (,+Af) , 

x n (t + At) = x„(t) + v„(t + At)At, 



(2a) 

(2b) 
(2c) 



where the function "Angle" is defined in such a way that for every vector u = ue' e we have Angle [5] = 9. In addition, <^„(f) 
is a random variable uniformly distributed in the interval and the noise intensity tj is a positive parameter. Note that 

for r\ = the dynamics are fully deterministic and the system reaches a perfectly ordered state, whereas for r\ > 1 the particles 
undergo purely Brownian motion. 

Analogously, the interaction rule with extrinsic noise is given by 



e n (t+At) 

v n (t + At) 
x n (t + At) 



= Angle 



y„(f) + T?^" w 



Voe ^(*+A0 

x n (t) + v n (t + At)At 



(3a) 

(3b) 
(3c) 



where £, n (t) and tj have the same meaning as before. Note that in Eq. d2al >. the angle of the average velocity V„(t) is computed 
and then a scalar noise TJ^„(f) is added to this angle. On the other hand, in Eq. (f3ab a random vector Tje'wW is added directly 
to the average velocity V„{t), and afterwards the angle of the resultant vector is computed. 14211 Intuitively, we can interpret the 
average velocity V n (t) as the signal that the n-th particle receives from its neighborhood, and the Angle function as the particle's 
decision-making mechanism. Therefore, in Eq. d2al l the particle clearly receives the signal from its neighborhood, computes the 
angle of motion, but then decides to move in a different direction. In contrast, in Eq. d3aT > the particle receives a noisy signal 
from the neighborhood, but once the signal has been received, the particle computes the perceived angle and follows it. 
At any time, the amount of order in the system is measured by the instantaneous order parameter defined as 



Wr,(t) 



1 



N 
n=l 



(4) 



where we have used the subscript 77 to explicitly indicate that the order in the system depends on the noise intensity. After a 
transient time, the system reaches a steady state, as indicated in Fig. [2] In the stationary state, we can average y/^ (f ) over noise 
realizations, or equivalently, over time. The stationary order parameter y/(r/ ) is, thus, defined as 



yr(fj) = lim(y»,(f)) 



1 f T 



(5) 
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FIG. 3: Phase transition diagram for a system with (a) extrinsic noise and (b) intrinsic noise. All the other parameters are the same for both 
systems: N = 2 x 10 4 , L = 32, vqAi = 0.05, and ro = 0.4. Note that in the extrinsic noise case the discontinuous phase transition is evident. 



where (•) denotes the average over realizations, whereas the average over time is given by the integral. 

Fig. [3] shows typical phase transition curves for both types of noise, intrinsic and extrinsic. The discontinuous character of 
the phase transition for the extrinsic noise case is evident (Fig. [3^). As far as we can tell, there are no doubts that the phase 
transition is discontinuous in this case, and there is plenty of numerical and theoretical evidence supporting this finding, l2j,|7|,|8|]. 
In contrast, the phase transition shown in Fig. [3J5 for the case with intrinsic noise seems to be continuous. However, despite 
this apparent continuity, strong finite-size effects, reflected in the smoothness of the curve as iff — > 0, may be hiding a possible 
discontinuity, as has been claimed by Gregoire and Chate, [2]. For some as yet unknown reason, these finite-size effects would 
be much stronger for the intrinsic noise than for the extrinsic noise, masking one discontinuity but not the other. 

The numerical determination of the nature of a phase transition is a difficult task precisely because the finite-size effects can 
conceal any discontinuity either in the order parameter or in its derivatives. Several techniques to deal with phase transitions 
numerically have been developed lf9il — H 13fl . For the part icular case of the Vicsek model, and other closely related models, finite- 
size scaling analysis and dynamic scaling analys is jl4ll-|[l8tl . renormalization group lfl9ll - ll22tl . mean- field theory 0], Q-JH, 
stability analysis j26ll . and Binder cumulants 121, 1271 l28tl . have been used. However, none of these techniques have determined 
conclusively the nature of the phase transition. 
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TABLE I: Parameters in the Vicsek Model 



Another problem for determining the nature of the transition, in addition to finite-size effects, is the existence of too many 
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FIG. 4: Phase transition in the Vicsek model for different system sizes. In all the cases the following parameters were used: p = 4, rg = 1, and 
vqAi = 0.03. Each point in each curve is the average of the instantaneous order parameter {jf-q (f ) over 10 6 time steps. 

parameters. This makes the numerical exploration of the parameter space very difficult, and what is valid for one region may not 
be valid for a different region. For instance, for low densities and small particle speeds, the phase transition looks continuous 
for systems as large as the computer capacity permits one to analyze, II 18l 12911 . But for high densities, large system sizes and 
high speeds the phase transition appears to be discontinuous, HQ. Table Exists all the parameters in the model. Of particular 
importance is / = (v'oAf ) / ro, the step size relative to the interaction radius. When I is small {I < 0.5) we will say that the system 
is in the low-velocity regime. Otherwise, the system will be in the high-velocity regime. 

III. THE CONTINUOUS PHASE TRANSITION 

The first claim of the existence of a continuous phase transition in the self-propelled particles model was presented in Vicsek's 
original work [Q]], which was entitled "Novel type of phase transition in a system of self-driven particles." The adjective "novel" 
was used probably because this phase transition was not expected in view of the Mermin- Wagner theorem 113011 , which asserts 
that in a low-dimensional (ID and 2D) equilibrium system with continuous symmetry and isotropic local interactions, there 
cannot exist a long-range order-disorder phase transition. Specifically, this theorem was proven having in mind Heisenberg and 
XY models of ferromagnetism. Vicsek's model in 2D is similar to the XY model, but in the latter the particles are fixed to the 
nodes of a lattice, whereas in the former they are free to move throughout the system. It is the motion of the particles that gives 
rise to effective long-range spatial interactions among the particles in the system. These interactions are responsible for the onset 
of long-range order, producing the phase transition. We will come back to this point later, (see Sec lIIIBl i. 

Fig.|4]shows the phase transition for different system sizes, keeping the density constant, p = 4. It should be emphasized that 
the curves shown in this figure correspond to a low-velocity regime in which the effective step-size / = v'oAf /ro is small: / = 0.03. 
This regime may be biologically plausible for systems such as flocks of birds, because birds can detect other birds that are far 
away compared to the distance traveled by each of them per unit of time.j43!] In Fig.|4j which is similar to the one originally 
published in Ref. [ 1 ] but with more particles, one can observe that the phase transition remains continuous even for very large 
systems, up to N = 131044. It is important to analyze very large systems because, as noted by Gregoire and Chate 121 1271 12811 . 
phase transitions occur in the thermodynamic limit defined as N — > °°, L — > °°, p = const. Using systems with different number 
of particles, up to N = 10 5 , Vicsek and his collaborators found numerically that the order parameter iff goes continuously to zero 
when the noise amplitude rj approaches a critical value r\ c from below as 

V~{ric-r]) P , (6) 

with a non-trivial critical exponent j3 = 0.42 ±0.03, iflil [l5ll . The critical value rj c of the noise at which the phase transition 
seems to occur depends on the system size L and on the density p . In this situation, finite-size scaling analysis can be done to 
check whether Eq. (O can be expected to be valid in the thermodynamic limit, and to determine if a universal function y/(rj /r\ c ) 
exists such that the relationship y/(rj) = \jr(rj / rj c ) holds for any size L and density p. 
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FIG. 5: Scaling of the phase transition with the system size, [panels (a), (b) and (c)], and with the density p, [panels (d), (e) and (f)]. In the 
three top panels the density was kept constant, p = 2, and the number of particles was increased. In the three bottom panels the system size 
was kept constant, L = 100, and the density varied, (p — 4 squares, p = 2 diamonds, p = 0.5 pluses). In all the cases the system was in the 
low-velocity regime / = 0.1. The figures were taken from Ref. 11411 . In that work the authors used rj £ [0, 2jt] and E, e [—0.5,0.5], whereas 
here we use rj 6 [0, 1] and | e [-7C, tc]. In (b) Arj = r\ c — r\. 



A. Finite-size scaling analysis 

In Refs. lfl4l[l5ll . Vicsek and his group presented a finite-size scaling analysis showing that a universal function \j/(rj / rj c ) does 
indeed exist in the low-velocity regime. Fig.|5]shows how the phase transition scales with the system size L and with the density 
p. From Figs.[5p and[5£ it is clear that the order parameter computed for different system sizes and different densities, plotted as 
a function of rj /r\ c , falls onto the same universal curve i//(tj / r\ c ). Additionally, the critical value of the noise rj c scales with the 
system size and the density respectively as rj c w L -1 / 2 , and r\ c ~ p K , with K ss 0.45. A similar scaling behavior was found for 
systems in ID and 3D lfl6l[l7ll . It should be noted that the scaling behavior reported in Fig.|5]is consistent with of second-order 
phase transitions. 

More recently, Baglieto and Albano performed a more thorough analysis IU8I1 . finding consistent evidence of a second order 
phase transition in the low-velocity regime (with / = 0.1) by means of two independent approaches. The first approach is 
finite-size scaling analysis, in which they measure stationary properties of the system such as the order parameter y/(T]) and 
the susceptibility % = (\j/ 2 ) — (y/) 2 . They assume that the following scaling relations, which are well established in equilibrium 
systems, hold for the Vicsek model close to the phase transition: 

y/(T7,L) « L-^>((7 7c -77)L 1 ' ,v ), (7) 

X (n,L) « L^((TJ c -rj)L 1/v ). (8) 

where y/ and % are the (universal) scaling functions. 

The second approach consists in performing short-time dynamic simulations and studying the results on the basis of dynamic- 
scaling theory, 13211 . The main idea is to start out the dynamics from fully ordered (V?j (?) = 1) and fully disordered (\f/^ (t ) w 0) 
configurations, and analyze how the system relaxes in time to the steady state. Close to the phase transition, the relaxation from 
an ordered state is characterized by the scaling relations l32tl 

y n {t,L) a b-Pl v \jf{b-Hy{r} c -r}),b- x L), (9) 
Xr,(t,L) w b^ v x(b- z ty(r, c ~r,),b- l L). (10) 

where z is the dynamic exponent, b is an arbitrary scale factor, and iff and % are the universal dynamic scaling functions. 




FIG. 6: (a) Schematic representation of the long-range spatial correlations developed in time in the Vicsek model due to the motion of the 
particles, (b) A network model in which long-range interactions have been explicitly introduced via a small-world topology. 



Additionally, it is known for equilibrium critical phenomena that the so-called hyperscaling relationship 

rfv-2j3 = y, (11) 

must be satisfied i33ll . Baglieto and Albano showed numerically the validity of Eqs.(f7l>-(fTTb for systems of different sizes 
and densities, obtaining excellent agreement between the numerical results and the expected scaling relationships. Thus, when 
finite-size effects are taken into account in the Vicsek model with intrinsic noise and in the low-velocity regime, two independent 
methods, one static and the other dynamic, show consistency with a second-order phase transition. 



B. Mean-field theory and network models 

As mentioned above, the XY model is similar to the Vicsek model in that neighboring particles tend to align their "spins" (or 
velocities) in the same direction. However, it is well known that in the XY model, in which the particles do not move but are 
fixed to the nodes of a lattice, there is no long-range order-disorder phase transition [33]. Thus, it is the motion of the particles in 
Vicsek's model that enables the system to undergo a phase transition. What this motion actually does is to generate throughout 
time long-range spatial correlations between the particles. This is illustrated in Fig. [6^, in which two particles that at time ?o do 
not interact, move and come together within the interaction radius at a later time t\. In Refs. ITJ and H23I1 a network model that 
attempts to capture this behavior is proposed. In this model the particles are placed on the nodes of a small-world network [34], 
in which each particle has K connections. The small-world property consists in that a fraction p of these connections are chosen 
randomly whereas the other fraction 1 — p are nearest-neighbor connections (see Fig. U). As in Vicsek's and the XY models, 
each particle is characterized by a 2D vector (or "spin"), v„(t) — voe' e "^\ of constant magnitude and whose direction changes in 
time according to Eqs. (O and ([TJ. However, in the network model the sum appearing in Eq. ([1) is replace by the sum over the 
K nodes connected to v„. 

The network model described above differs from the Vicsek model in that the particles do not move, but are fixed to the 
nodes of a network. However, the long-range spatial correlations created by the motion of the particles in Vicsek's model are 
mimicked in the network by the small-world topology, namely, by the pKN random connections scattered throughout the system. 
Although the analogy is not exact, the network model can be considered as the mean-field version of Vicsek's model, particularly 
when p = 1, i.e. when the K connections of each particle come randomly from anywhere in the system. In this case all spatial 
correlations are lost and the mean-field theory is exact. In Ref. B23I1 it is shown analytically that the network model with intrinsic 
noise and p = 1 undergoes a second-order phase transition. In contrast, in Refs. 0, l24ll it is proved also analytically, and for 
p = 1, that the network model with extrinsic noise exhibits a first-order phase transition. In Ref. [7] it was shown numerically 
that the network model with p = 1 corresponds to the Vicsek model in the limiting case of very large particle speeds (l^S> 1). This 
is because of the fact that for very large particle speeds, two neighboring particles that at time to are within the same interaction 
vicinity, because of the noise in the direction of motion of each particle, at the next time step fo + At, they will end up far away 
from each other and interact with other particles. Thus, in the limit vq — > °° of the Vicsek model, all spatial correlations are 
destroyed in one single time step, which is exactly the condition for the mean-field theory to be exact. 
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FIG. 7: Phase transition in the network model with (a) intrinsic noise and (b) extrinsic noise. Note that in the first case the phase transition is 
continuous whereas in the second case it is discontinuous with a region of hysteresis. The network models can be considered as the mean-field 
description of the Vicsek model, especially for p = 1. However, the same behavior is observed for any non-zero value of p, even for very small 
ones. For p = the small world network model is analogous to the XY model and no long-range order phase transition exists. 



Numerical simulations show that for any p S (0, 1], the phase transitions with intrinsic noise is continuous, whereas it is 
discontinuous with extrinsic noise (see Fig.|7]i. The results obtained with the network model indicate that the way in which the 
noise is introduced into the system can indeed change the nature of the phase transition. Therefore, the fact that the extrinsic 
noise produces a first-order transition cannot be used as supporting evidence for the continuous or discontinuous character of the 
transition with intrinsic noise. [44] 

The mean-field approach conveyed by the network model has been criticized on the basis that it neglects the coupling between 
local density and global order produced by the motion of the particles, ll27ll . This coupling can generate spacial structures such 
as clusters of particles I331 . whereas in the mean-field theory one assumes that the system is spatially homogeneous. While this 
is true, is not clear how the coupling between order and density in the Vicsek model changes the nature of the phase transition. 



C. Hydrodynamic models 

For relatively large densities, one can imagine a flock as a continuum fluid in which neighboring fluid cells tend to align their 
momenta. This was the approach taken by Toner and Tu, lfl9tl - l22ll . who constructed a phenomenological hydrodynamic equation 
which has the same symmetries as the Vicsek model. Specifically, this equation includes isotropic short-range interactions, a 
self-propelled term, a stochastic driving force, and does not conserve momentum. Together with mass conservation, the system 
of hydrodynamic equations studied by Toner and Tu is 

<9 r v + ^i(y- V)v + ^2(V-v)y + A 3 V(|v| 2 ) = av - J3 |vfv - VP + D B V(V ■ v) +D T V 2 v + D 2 (v ■ V) 2 v + f, (12) 

P(P) = £>„(p-po)", (13) 

n=\ 

d t p+V-(vp) = 0, (14) 

where v is the velocity field, p the density, P the pressure, rj„ are the coefficients in the pressure expansion, and / is a random 
Gaussian white noise which is equivalent to the noise term in the Vicsek model. Through a renormalization analysis of the 
above equations, Toner and Tu found scaling laws and relationships between different scaling exponents. Although the authors 
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do not explicitly mention the nature of the phase transition, the characterization they present in terms of scaling laws and scaling 
exponents is, again, not inconsistent with a continuous phase transition in 2D. 

Later, Bertin et al. derived an analogous hydrodynamic equation from microscopic principles incorporating only binary 
"collisions" of self-propelled particles with intrinsic noise into a Boltzmann equation, j26ll . One of the advantages of this 
approach is that all the phenomenological coefficients appearing in Eq. ( fT2l . (the A's, 0£, j3, etc.), can be explicitly expressed in 
terms of microscopic quantities. The authors then show that the hydrodynamic equation has a homogeneous non-zero stationary 
solution that appears continuously as the density and the noise are varied. However, linear stability analysis shows that for this 
model, this homogeneous solution is unstable and therefore, the authors conclude, more complicated structures should appear in 
the system (such as bands or stripes, see Sec. Ml. 

In any case, the evidence collected so far with hydrodynamic models of swarming, formulated along the same lines as the 
Vicsek model, while is not inconsistent with a second order phase transition, it has not been able either to shed light on the nature 
of the transition. 



D. Models based on alignment forces 

In a different type of swarming models the interactions between the particles are modeled through forces that tend to align 
the velocities of the interacting particles. Newton equations of motion are then solved for the entire system. Along these lines, 
Peruani et al. analyzed a system of over-damped particles whose dynamics is determined by the equations of motion Il25ll 

dx m (t) _ je m (t) 



dt 



= v e lw > (15a) 



djjt) = _ Y du m (* e n ) 

dt a a m 

where 7 is a relaxation constant, r\ m (t) is a Gaussian white noise (similar to the intrinsic noise in Vicsek's model), and the 
potential interaction energy U m (x n , Q n ) is given by 

U m {x n , 0„) = cos (0m - Qn)- 

\x m -x„\<r 

where ro is the radius of the interaction vicinity. In the limiting case 7 — > °° of very fast angular relaxation, the above model is 
completely equivalent to the Vicsek model, B25I1 . 

Analogously, inspired by the migration of tissue cells (keratocytes), Szabd et al. formulated a model of over-damped swarming 
particles for which the dynamic equations of motion are l36ll 

dx m (t) N - 

v m {t) = — 7 — =von m {t) + liYF(x m ,x„) (16a) 
dt , 

n=l 



dejt) 1 

= — arcsin 

dt T 



n m (t) x 



(01 



•TJ m (0 (16b) 



where m is the angle of the unit vector n m , and rj„,(f) is again a Gaussian white noise. The interaction force F(x m ,x„) depends 
piece-wise linearly on the distance d nm = \x m —X n \. It is repulsive if d mn < Rq, attractive if ro < d mn < R\, and vanishes if 
Ri < d mn . The explicit form of the force is 

f F r ±^ if d mn < Ro 

F{x m ,x n ) = e mn x \ F a ^f^ if R < d m „ < Ri 

I if Ri < dmn 

where e m „ is the unit vector pointing from x„ to x,„, and F r and F a are the magnitudes of maximum repulsion and maximum 
attraction, respectively. 

In the two models described above, the forces between neighboring particles tend to align their velocities if the particles are 
withing the specified region, and in both models the particles are subjected to intrinsic noise. Interestingly, numerical simulations 
show that the states of collective order in these two models appear continuously through a second order phase transition, as Fig. [8] 
illustrates. 
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FIG. 8: Phase transitions in swarming models based on alignment forces with intrinsic noise, (a) Order parameter as a function of the noise 
intensity for the model given in Eqs. d 1 5 b - The dashed curve is the mean-field prediction whereas the symbols correspond to numerical 
simulations, (b) Phase transition for the model given in Eqs. d!6t . The data were obtained from numerical simulations. In both cases the phase 
transition is continuous. The figures (a) and (b) were taken (with small modifications to make the labels clearer) from Refs. 1231 and I3ql . 
respectively. 



IV. THE DISCONTINUOUS PHASE TRANSITION 



The first evidence for the discontinuous character of the phase transition in the Vicsek model with intrinsic noise appeared in 
Ref. (2], and it was later reinforced in Ref. lf28ll . In that work the authors studied very large systems in the high-velocity regime 
I > 0.5, 14511 . Fig. [9] summarizes this evidence for a large system, with N = 131072 particles, L = 256, and / = 0.5. This system 
is similar to one of the systems used by Gregoire and Chate in Ref. [2] to demonstrate the discontinuity of the phase transition 
with intrinsic noise. It is clear from Fig. [9^ that the order parameter undergoes an apparently discontinuous phase transition 
around a critical value r\ c = 0.478. Such a discontinuity is confirmed by the Binder cumulant G(r\), defined as 

CM = 1-5^- ("> 

The Binder cumulant, which is closely related to the Kurtosis, can intuitively be interpreted as a measure of how much a 
probability distribution function (PDF) differs from a Gaussian, particularly concerning the "peakedness" of the distribution. It 
was introduced by Binder as a means to determine numerically the nature of a phase transition, |[T2il39ll . The idea behind this 
measure is that in a first-order phase transition, typically two different metastable phases coexist. Throughout time, the system 
transits between these two phases and consequently the time series of the order parameter shows a bistable behavior, as it is 
illustrated in Fig. [9};. The PDF P(y/jj) of this time series is then bimodal, reflecting the coexistence of the two different phases 
and the transit of the system between them (see Fig.[9jl). Because of this bimodality of which occurs only close to the 

phase transition, the Binder cumulant has a sharp minimum, as it is actually observed in Fig. [9^. On the other hand, far from the 
phase transition, either in the ordered regime or in the disordered one, P(y/jj) is unimodal and rather similar to a Gaussian, for 
which the Binder cumulant is constant. It is only close to a discontinuous phase transition that P{^rf) becomes bimodal, thus 
deviating from a Gaussian, and the Binder cumulant G(rj) exhibits a sharp minimum. This is the strongest evidence presented 
by Gregoire and Chate in favor of a discontinuous phase transition in the Vicsek model with intrinsic noise, lHIH. They also 
analyzed variations of this model, for instance by introducing cohesion or substituting the intrinsic noise by extrinsic noise. 
However, the mean field network models presented in Sec. IIIIBl show that the results obtained for one type of noise cannot be 
extrapolated to a different type of noise. Thus, in this section we will focus on the intrinsic noise case. 

It should be emphasized that the discontinuous behavior displayed in Fig. [9] is observed only in the high velocity regime. 
Indeed, it is apparent from Fig. [4] that in the low-velocity regime I = 0.03, the order parameter y/(rj) does not exhibit any 
discontinuity even for a large system with N = 131044 and L = 181, (comparable in size to the one used in Fig. [9]). It could be 
possible that the phase transition changes from continuous to discontinuous across the parameter space. However, the authors of 
Refs. 10] and j28fl claim that this does not happen. They state that all phase transitions occur in the thermodynamic limit, and if 
we were able to explore very large systems, then the phase transition in the Vicsek model with intrinsic noise would turn out to 
be discontinuous for any values of the parameters, M27I1 . 

Another important aspect of the discontinuous phase transition shown in Fig. [9^ is the absence of hysteresis. Indeed, our 
numerical simulations show that for this system it is irrelevant to start out the dynamics from ordered or from disordered states, 
because in either case the curve displayed in Fig. [9^ is obtained. This is in marked contrast with the hysteresis commonly 
observed in usual first-order phase transitions, as the one depicted in Fig. [7j), which corresponds to the extrinsic noise dynamics 
and therefore is undoubtedly first order and in accord with the mean-field analytical results (see Sec. lIIlBT i. In the next section 
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FIG. 9: Evidence of the discontinuous character of the phase transition in the Vicsek model with intrinsic noise, for as system with N = 131072 
particles, L = 256, and / = 0.5. (a) Order parameter \jf as a function of the noise intensity tj. In this case the order parameter decays to zero 
discontinuously around the critical value r\ c w 0.478. (b) Binder cumulant G as a function of TJ. Note the sharp minimum of the Binder 
cumulant close to the phase transition, (c) Time series of the instantaneous order parameter iff^ (t) for tj = 0.475 (very close to the phase 
transition). This time series exhibits bistability, which is characteristic of phase coexistence in a first-order phase transition, (d) Probability 
distribution function f (1/67) of the instantaneous order parameter \jfrj (?). This PDF is bimodal because of the bistability of the time series. In 
(a) and (b) each data point is the average of y/jj (?) over 2 x 10 7 time steps. 



we present evidence indicating that the discontinuity of the phase transition observed in Fig. |9]may in fact be an artifact of the 
periodic boundary conditions, as first proposed by Nagy et al, ll29ll . 



V. THE EFFECT OF THE BOUNDARY CONDITIONS 



In a recent paper B29H . Nagy, Daruka and Vicsek pointed out the formation of traveling bands in the high velocity regime. In 
Fig. [10] we show these bands for three different values of the noise intensity tj and for the same system used to compute the 
phase transition in Fig. [9] Note that the three values of tj in Fig. [Tol (tj = 0.430, TJ = 0.440 and tj = 0.450), correspond to 
the ordered regime, quite far from the transition, which occurs around tj 6 w 0.478. Interestingly, the bands always wrap around 
tightly in a periodic direction, which for simple shapes is usually either perpendicular or parallel to a boundary of the periodic 
box, or along a diagonal. In the ordered regime, the bands are very stable. Once they form, they may travel for extremely long 
times in the same direction (as shown in Fig.QJJ the motion of the bands, indicated by the arrow, is perpendicular to themselves). 
In Ref. M29I1 the authors also use hexagonal boxes with periodic boundary conditions, and the bands again always form parallel 
to the boundaries of the hexagon or to its diagonals. 

When the value of tj approaches the critical value tj c , the bands become unstable: They appear and disappear intermittently 
throughout the temporal evolution of the system. Fig. [TJJshows the time series of the order parameter for the same system used 
in Fig. [9] at a noise intensity tj = 0.475, i.e. very close to the discontinuity. The time series shows that the order parameter 
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FIG. 10: Traveling bands in the Vicsek model with intrinsic noise. The top panel shows the time series of the order parameter for three different 
values of the noise intensity T7. The images at the bottom are snapshots of the system corresponding to each of the three values of 77 used. The 
arrows indicate the direction of motion of the bands. 



"jumps" between two rather different values. These jumps give rise to the bimodality of the order parameter, which, according 
to the interpretation presented in the previous section, indicates phase coexistence. A close inspection of the spatial structure of 
the system reveals that such bistability is produced by the creation and destruction of the traveling bands, as the two snapshots in 
Fig.[TTJindicate. Thus, as time passes, the system transits from spatially homogeneous states to states characterized by traveling 
bands. The authors of Ref. 112811 interpreted these bands as one of the two metastable phases, which coexists with a more or 
less homogeneous background of diffusive particles (the other metastable phase). However, Vicsek and his group retorted that 
the formation of traveling bands is not an intrinsic property of the system, but, rather, an artifact of the periodic boundary 
conditions. i46ll 

It is a very common resource in numerical simulations of physical systems to use periodic boundary conditions as way to 
represent an infinite homogeneous space and mitigate boundary effects. In Vicsek's model the particles move in a box of sides 
L. If a given particle exits the box through one of the sides, it is re-injected into the system through the opposite side at a 
specular position (see Fig. [T2"h). Under these circumstances, the real geometry of the box is that of a torus, and the motion 
of the particles takes place on its surface. The toroidal topology induced by the periodic boundary conditions clearly favors 
the formation of bands, which in fact are closed stable rings moving on the surface of the torus, as is illustrated in Fig. [T2b . 
However, there is no reason whatsoever to re-inject the particles at specular positions of the box. Rather, it is possible to re-inject 
the particles at positions that are shifted a distance AL with respect to the specular position, as is done in Fig.[T2b. By doing so, 
the toroidal topology is broken and the formation of bands (closed rings) is not longer stabilized by (toroidal) periodic boundary 
conditions. The idea is schematically illustrated in Fig. [T2t . although the real situation is more complicated because the shift 
AL is implemented in the four sides of the box (see the Appendix). This shift, which we will refer to as the boundary shift, does 
not affect the homogeneity of the system, as all the particles are subjected to it when they hit the boundary. Interestingly, when 
the shift is implemented, the discontinuity of the phase transition reported in Refs. |01 and l28tl . and reproduced here in Fig. [9] 
seems to disappear. 

Fig. [T"3h shows the phase transition for the same system as the one used in Fig. [9] but this time comparing the two cases: 
With fully periodic boundary conditions, i.e. AL = (black curve with circles), and with a non-zero boundary shift AL = 0.2L 



13 




le+06 2e+06 3e+06 4e+06 5e+06 6e+06 7e+06 

T 




B 




W^mmm 



FIG. 11: Bistable time series of the order parameter at a noise intensity rj = 0.475 for the same system as in Fig. [9] This value is very close to 
the critical value rj c = 0.478. The bottom figures are snapshots of the system at the times indicated by the points A and B in the time series. 
The apparent bistability of the order parameter indicates the formation and destruction of the bands. 



(red curve with squares). In Fig. [T3b the corresponding Binder cumulants are reported. We choose AL = 0.2L because this 
is approximately the width of the bands close to the transition (see snapshot A in Fig. fTTT i. It is clear from Fig. [13] that the 
introduction of the boundary shift eliminates the discontinuity of the phase transition, making it continuous. Note that the sharp 
peak in the Binder cumulant also disappears. It could be argued, then, that if the phase transition was intrinsically discontinuous, 
the introduction of the boundary shift should not change its character. This is indeed what happens in the Vicsek model with 
extrinsic noise with a boundary shift AL = 0.2L, as illustrated in Fig. [14] In such case, of course, the introduction of the 
boundary shift does not change the appearance of the phase transition. The above results give support to the claim that the 
apparent discontinuous character of the phase transition in the intrinsic noise case may indeed be an artifact of the topology 
induced by the periodic boundary conditions, as Vicsek and his colleagues had already suggested, [29]. 

A more dramatic effect of the periodic boundary conditions on the phase transition can be observed in a simple variation of 
the Vicsek model, for which it is not necessary to analyze very large systems. In this variation, the particles have a "vision 
angle". In other words, the interaction region of each particle is no longer a full circle surrounding it, but only a circular cone 
of amplitude 2a centered at the particle's velocity, as illustrated in Fig. [T5k All the other aspects of the interaction are as in 
Vicsek's model, only the form of the interaction vicinity is changed. In what follows we consider such a system with intrinsic 
noise and N = 1600, L = 20, 1 = 0.03 (low-velocity regime), and a = OAk. For periodic boundary conditions and some values 
of the intrinsic noise, after a transient time the particles form very dense bands, as shown in Fig. [15b and the corresponding 
snapshots. Actually, as time goes by, the system transits back and forth, from configurations like the one depicted in snapshot 
A, to configurations with bands, as the one shown in snapshot B. Those bands produce two discontinuities in the curve of the 
order parameter, which is plotted in Fig.fTBtl (black curve with circles). Interestingly, these discontinuities do not occur close to 
the phase transition (which in this case seems to happen around r] c m 0.6), but well into the ordered phase. On the other hand, 
Fig. [TBb shows the time series of the order parameter when a boundary shift AL = 0. 1L is introduced. Note that in this case the 
order parameter does not reflect the transit of the system between any two phases. This is because with the boundary shift the 
bands no longer form and the two discontinuities in the phase transition curve disappear, as is apparent from Fig.[l5tl (red curve 
with squares). 

This variation of the Vicsek model clearly demonstrates that the periodic boundary conditions can have a dramatic effect 
on the system dynamics, producing artificial discontinuities in the order parameter that would not necessarily be present in an 
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FIG. 12: (a) In the original Vicsek model one implements periodic boundary conditions: If a particle exits through one side of the box, it 
re-enters through the opposite side at a specular position, (b) The periodic boundary conditions induce a toroidal topology that favors the 
formation of bands (closed rings) moving on the torus' surface, (c) Shift AL in the boundary conditions. The shift is implemented in the four 
sides of the box, which breaks the toroidal symmetry of the dynamics and makes the closed ring structures unstable, as it is illustrated in (d). 

infinite system with no boundaries. 

VI. SUMMARY 

In this work we have presented much of the evidence offered in support of two distinct points of view, at odds with each other, 
about the nature of the phase transition in the Vicsek model of self-propelled particles with intrinsic noise. On the one hand, there 
is theoretical and numerical evidence indicating that such a phase transition is continuous. This evidence comes from several 
fronts: From the numerical side, finite-size scaling analysis and dynamic scaling analysis performed for systems operating in the 
low-velocity regime, have yielded results consistent with a continuous phase transition. These numerical techniques show scaling 
behavior of the order parameter with the size of the system and with the density, for systems as large as the computer capacity 
allows one to analyze. Additionally, the hyperscaling relationship between the critical exponents has been confirmed within 
numerical accuracy. This scaling behavior is consistent with what is known for second order phase transitions in equilibrium 
systems. From the theoretical side, analytic calculations unambiguously show that the mean-field version of the Vicsek model 
with intrinsic noise undergoes a second-order phase transition. Furthermore, swarming models based on alignment forces, (rather 
than on kinematic rules), also exhibit continuous phase transitions when the intrinsic noise is used. 

On the other hand, there is direct numerical evidence suggesting that the phase transition is discontinuous in very large 
systems operating in the high-velocity regime. In addition to a clear discontinuity in the phase transition curve for these large 
systems, this point of view is supported by an apparent bistability of the order parameter close to the phase transition, and the 
corresponding bimodality of its PDF signaled by a "dip" of the Binder cumulant. Such a scenario is characteristic of phase 
coexistence in first-order transitions. 

However, numerical computations performed for systems comparable in size to those used to demonstrate that the phase 
transition is discontinuous, suggest that the apparent discontinuity is an artifact of the toroidal topology induced by the periodic 
boundary conditions. For, when the boundary conditions are "shifted" in order to break the periodicity, and hence the toroidal 
symmetry, the discontinuity disappears. This effect is enhanced when the particles have a limited vision angle. 

Determining the nature of the phase transition in the Vicsek model with intrinsic noise is important, among other things, 
for understanding the physical principles through which collective order emerges in swarms and flocks, and more generally, in 
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FIG. 13: Effect of the boundary shift AL on the phase transition in the Vicsek model with intrinsic noise, (a) Order parameter as a function 
of the noise intensity for a system with no shift (black curve with circles), and for a system with a boundary shift AL = 0.2L (red curve with 
squares). Note that the apparent discontinuity of the phase transition for the system with periodic boundary conditions disappears when the 
shift is introduced. " N " ' ' ' '' ' '' ' '' ' ' ' ''' ' '' ' '' the system 
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FIG. 14: Phase transition in the Vicsek model with extrinsic noise for systems with N = 80000, L = 64, and / = 0.5. The black circles 
correspond to fully periodic boundary conditions, i.e. AL = 0, whereas the squares correspond to a boundary shift AL = 0.2L. It is clear that 
in this case the boundary shift does not change the phase transition, which is undoubtedly discontinuous. 



systems out of equilibrium. Despite the efforts of many groups during the last decade, nobody has said the last word about the 
onset of collective order in swarming systems. This remains an open problem. 
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FIG. 15: Phase transition in a variation of the Vicsek model with N = 1600, L = 20 and / = 0.03. (a) The interaction vicinity is a circular 
cone of radius ro = 1 and angular amplitude 2a, centered at the particle's velocity, (b) Time series of the order parameter for a = OAk and 
77 = 0.44, with fully periodic boundary conditions (AL = 0). Note that after a transient time, the particles aggregate in a very dense band 
parallel to the boundary. The snapshots A and B correspond to the points signaled with the arrows in the time series, (c) The same as in the 
previous panel but now with a boundary shift AL = 0. 1L. (d) Phase transition diagram with (red curve) and without (black curve) boundary 
shift. Note the discontinuities produced by the periodic boundary conditions. 

APPENDIX A: BOUNDARY SHIFT 

In this appendix we explain how the boundary shift AL was implemented in our numerical simulations. Let us start with a 
system with periodic boundary conditions, as in Fig. [12b. Using periodic boundary conditions is equivalent to making infinite 
copies of the original system, and covering homogeneously all the space with these copies. This is illustrated in Fig.[T6k where 
a central box has been "cloned" several times. A particle exiting the central box just re-enters into it at a position equivalent to 
the one it would have in the adjacent box, and so on. 

On the other hand, Fig. [16b shows how the boundary shift AL is implemented. The copies of the central box have been 
shifted around it the same positive distance AL in both the horizontal and the vertical directions. In this case, the space is not 
cover homogeneously because the shifted boxes overlap, as it is clearly seen in Fig. [16b . However, this overlap is important 
only for particles exiting the central box from the corners (as the one signaled with an arrow), because such a particle does not 
"know" what box to go to. This overlapping mixes the particles exiting from the corners in a way which is different from the 
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FIG. 16: (a) A system with periodic boundary conditions. This is equivalent to making infinite copies of the central box and covering 
homogeneously the space with them, (b) A system with boundary shift. Note that the shifted boxes overlap at the corners of the central box. 



natural motion of the particles located at the bulk of the system. Thus, the boundary shift AL introduces inhomogeneities in 
the "tiling" of the space. These inhomogeneities occur in regions of size VqA/ around the corners of the central box. Therefore, 
if (vrjAf)/L <C 1, as it has been the case in all the simulations carried out in this work, the inhomogeneities generated by the 
boundary shift should be negligible. In our simulations, when a particle exits the central box from one of the corners, it randomly 
chooses one of the two overlapping boxes to be re-injected into it. 
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